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SECOND  ORDER  APPROXIMATION  METHOD  FOR  TRANSONIC  FLOWS  OVER 
SWEPT  WINGS 


Shen  Keyang 

Shanghai  Aviation  Industry  Corporation 
Abstract 

This  paper  uses  the  TSDH  equation  to  calculate  the  transonic 
flows  over  three-dimensional  swept  wings.  It  considers  the 
leading  edge  boundary  conditions  and  leading  edge  velocity 
potential  equation  suitable  for  the  blunt  leading  edge  of  wings. 
We  use  the  expanded  form  of  the  Jameson  scheme  in  the  unequal 
step  length  grid  for'  the  disretization  of  the  TSDH  equation  into 
a  set  of  finite  difference  equations.  Afterwards,  we  arranged 
a  rarefied  grid  in  the  calculation  space,  then  arranged  a 
dense  grid  near  the  wing  and  carried  out  alternate  iteration  of 
the  rarefied  and  dense  meshes  so  as  to  accelerate  convergence 
and  raise  calculation  precision.  Calculations  of  the  super¬ 
critical  without  shock  wave  and  with  shock  wave  conditions  of 
the  ONERA  M6  wing  show  that  there  is  good  agreement  between  the 
TSDH  solution  and  the  FVP  solution  and  wind  tunnel  tests. 

Symbols 

Q  :  local  velocity  of  sound 
E  :  iteration  error 
M  :  Mach  number 
n  :  iteration  number 
a  :  attack  angle 
Y  s  specific  heat  ratio  of  air 
6  :  mean  relative  thickness  of  wing 
£  :  small  disturbance  parameters,  £  =  <5  2//^/Mcc 

/i  :  sweepback  angle 
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Lower  Symbols 


«>  :  free  flow 
LE  :  leading  edge 
TE  :  trailing  edge 


Upper  Symbols 


(n) 

x,y,z 

Q 

y  (x) 
■*w 


M' 

qp 


u 

1 

(n-1) 


:  the  n  number  iterations 

:  right  angle  coordinate  system,  x  axis  along  wing 
chord  direction,  the  y  axis  goes  up 

:  flow  velocity 

:  wing  configuration  coordinate 

:  switch  parameters 
:  full  velocity  potential 
:  disturbance  velocity  potential 
:  damping  coefficient 

:  low  relaxation  coefficient 

:  upper  wing  surface 
:  low  wing  surface. 

:  the  (n-1)  number  iterations 


I.  Introduction 


Because  of  the  simplicity  of  the  calculations  of  the  classi¬ 
cal  transonic  small  disturbance  theory  (TSD  equation) ,  it  has 
continuously  been  given  attention.  However,  when  compared  with 
the  full  velocity  potential  theory  (FVP  equation) ,  it  has  the 
following  two  major  drawbacks:  (1)  because  a  longitudinal 
disturbance  velocity  component  is  not  designed  into  the  wing 
surface's  boundary  conditions,  infinitely  great  singularity  is 
produced  in  the  area  of  the  blunt  leading  edge;  (2)  based  on 
the  conditions  of  the  existence  of  shock  waves  on  the  infinite 
swept  wing  [1] ,  the  TSD  equation  is  only  suitable  for  situations 
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wherein  the  wing’s  sweepback  angle  is  relatively  small  (see 
Appendix  A) .  People  hope  to  attain  a  calculation  method  which 
will  on  the  one  hand  maintain  the  simplicity  of  the  TSD  theory 
and  can  also  overcome  the  two  above  mentioned  drawbacks.  Ballhaus 
and  Bailey  [1]  used  the  MSD  equation  to  basically  overcome  the 
above  drawback  (2)  yet  it  does  not  consider  drawback  (1) .  In 
recent  years,  References  [2,3]  researched  the  second-order 
approximation  method  (TSDH  equation)  in  the  transonic  small 
disturbance  potential  flow  and  its  application  fob  the  flow  over 
airfoil  overcame  the  above  drawback  (1)  and  obtained  results 
very  close  to  the  FVP  solution.  The  aim  of  this  paper  is  to 
extend  this  method  to  the  transonic  flow  over  three-dimensional 
wings  and  thus  overcome  the  above  drawbacks  of  the  TSD  equation 
and  obtain  results  very  close  to  the  FVP  solution. 

In  order  to  be  able  to  be  suitable  for  calculation  of  wings 
with  complex  configurations  as  well  as  wings  having  a  fence, 
leading  edge  sawteeth  or  wingtip  winglet  etc.  aerodynamic 
equipment,  this  paper  does  not  carry  out  certain  coordinate 
conversion  of  the  TSDH  equation  but  directly  solves  its 
difference  equation  in  the  physical  space.  At  the  same  time, 
this  paper  uses  Boppe's  [4]  rarefied  and  dense  mesh  alternate 
iteration  method  to  raise  the  convergence  speed  and  reduce  the 
requirements  for  computer  storage.  Calculations  of  ONERA  M6 
wing  with  complex  shock  wave  form  show  that  the  method  in  this 
paper  is  effective. 

II.  Basic  Equations  and  their  Boundary  Conditions 

Based  on  the  research  of  References  [2,3],  the  second-order 
and  some  of  the  third-order  terms  of  the  transonic  small 
disturbance  velocity  potential  equation  (TSDH  equation)  retained 
in  the  three-dimensional  situations  can  be  expressed  as  follows: 
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In  the  formula,  M'  =1  and  when  /J*  =0,  formula  (1)  is  simplified 
into  the  TSD  equation.  The  wing  surface's  boundary  conditions 
(except  for  the  leading  edge)  can  be  expressed  as 

cp,(jc,  ±0,  *)-si  !(^r)u.rC09a  +‘4*'p«<  *.  tO,  z))-sinol  (2) 

In  the  formula,  £  =f  /£  ,  f  is  the  correction  coefficient  and 
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this  paper  uses  f  -M^  (when  ^  1) .  In  the  blunt  leading 
edge  x=x  (z)  area,  (t?y  /<?x)  •,  =  «»  and  its  leading  edge 
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boundary  condition  is 
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We  substitute  formula  (3)  into  the  FVP  equation 
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(1) 

The  obtainable  simplified  FVP  equation  of  the  wing's  blunt 
leading  edge  is 
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In  the  formula 
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Because  we  only  consider  the  longitudinal  flow,  there  are 
symmetric  conditions  on  the  wing's  central  line 


(7) 

The  sweepback  discontinuity  surface  conditions  and  distant  field 
boundary  conditions  of  the  wing  are  identical  to  those  in 
Reference  [5].  Wing  surface  pressure  coefficient  Cp  is  calculated 
with  the  Bernoulli  equation: 

--^=^-A/iC2E<P.+c*<P2+^(20E<Py+eI<p;+e1qpJ)3lT/,T  lj  (8) 

III.  The  Finite  Diffference  Form  and  Solution  Equation 

Based  on  the  discussion  in  Reference  [3]  regarding  numerical 
calculation  stability  and  artificial  viscosity  directionality, 
formula  (1)  only  has  which  should  based  on  the  changes  of 

the  equation  form  use  the  central  difference  or  upstream 
difference  form  and  each  of  the  rest  of  the  terms  then  always 
use  the  central  difference  form.  Therefore,  the  difference  form 
of  this  method  is  non-r'otational  and  the  calculations  are  simpler 
than  those  of  the  full  velocity  potential  equation.  We  use  the 
extended  form  of  the  Jameson  scheme  in  a  non-uniform  grid.  For 
example,  in  supersonic  points  (i,j,k),  we  can  change  formula  (1) 
into 
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In  the  formula,  we  take  CC  ,  ^0,  p=x .  . -x .  _ ,  q=x.-x.  ., 
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is  calculated  by  the  upstream 


difference  formula. 
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In  the  formula,  SP  ^  as  well  a's  ’■P  and  ^  are  all 

calculated  by  the  central  difference  formula.  After  solving 

!n^.  by  the  set  of  difference  equations,  we  can  use  the 
•  1  t  J  r  K 

following  formula  to  ob.tain  fnl  . 

1  f  J/K 

(10) 

In  the  formula,  we  take  tO  ^  1. 

P 

For  each  term  in  leading  edge  equation  (5) ,  <p„„,  ®,  and 

xx  xy 

t  zz  use  a  one-sided  difference  formula  to  satisfy  the  boundary 

conditions  of  formula  (3)  and  V  and  then  use  the 

yy  zz 

central  difference  formula.  Moreover,  only  when  /U^=l  can 
formula  (5)  participate  in  the  operation  and  combine  with 
formula  (1)  for  solution. 

When  solving  velocity  potential  (p  ,  we  arrange  two  grids 
in  the  calculating  space:  we  arrange  a  rarefied  mesh  in  the 
calculating  space  (for  example,  30x21x15,  at  this  time  the 
spanwise  direction  of  the  wing  surface  is  cut  into  several 


sections  and  each  section  has  5-7  chordwise  mesh  points)  and  at 
the  same  time  we  arrange  a  dense  mesh  near  the  wing  (the  density 
of  the  dense  mesh  is  fixed  in  view  of  requirements,  the  most 
rarefied  situation  is  35x25x15  and  at  this  time  the  spanwise 
direction  of  the  wing  surface  is  cut  into  more  than  ten  sections 
and  each  section  has  15-20  chordwise  mesh  points) .  Afterwards, 
we  carry  out  rarefied  and  dense  mesh  alternate  iteration.  The 
solution  process  is  carried  out  in  three  stages: 

The  first  stage:  we  give  a  velocity  potential  initial  field 

in  the  rarefied  mesh,  for  example,  we  take  .  =0.  Taking 

i  /  D  /  k 

/t^=0,  we  carry  out  iteration  solution  of  equations  (1)  and 
and  at  the  same  time  continuously  change  the  distant  field 
velocity  potential  with  the  increase  of  the  circulation. 

Because  the  grid  is  relatively  rar'efied,  after  a  relatively 
small  number  of  iterations,  it  can  reach  convergence  standard 

E=MAX  /  <p  10  2-10  3.  In  this  way,  we  establish 

the  approximate  velocity  potential  field  and  circulation 
distribution  in  the  entire  calculating  space. 

The  second  stage:  this  stage  is  carried  out  in  two  steps. 

The  first  step  is  finding  the  initial  field  of  the  dense  mesh 

velocity  potential  from  the  rarefied  mesh  velocity  potential 

interpolation.  For  the  velocity  potential  value  of  the  fixed 

dense  mesh's  outer  boundary,  we  carry  out  iteration  solution 

using  in  a  dense  mesh  for  formulas  (1)  and  (2)  so  as  to 
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attain  convergence  standard  E  ^  10  -10  .  The  second  step  is 

to  find  the  velocity  potential  value  on  the  wing  surface  in  the 
rarefied  mesh  from  the  dense  mesh's  velocity  potential  interpo¬ 
lation  found  in  the  first  step,  to  find  the  rarefied  mesh's 
distant  field  velocity  potential  value  from  the  dense  mesh's 
circulation  distribution  and  the  two  remain  unchanged  during  the 
iteration  process.  Therefore,  this  step  solves  the  rarefied 
mesh's  velocity  potential  under  Dirichlet  boundary  conditions. 


At  this  time,  we  still  use  A*=o  and  after  attaining  the 
convergence  standard  we  again  repeat  the  first  step.  After 
alternate  iteration  of  the  rarefied  and  dense  meshes  and  the 
wing  surface  pressure  coefficient  in  the  dense  mesh  reaches  a 
certain  convergence  standard,  we  can  conclude  the  operation  of 
the  second  step. 

The  third  step  is  carrying  out  densif ication  of  the  dense 
mesh  so  as  to  attain  the  number  of  mesh  points  sufficient  for 
the  needs  of  the  wing.  At  this  time,  the  velocity  potential 
value  of  the  fixed  dense  mesh's  outer  boundary  uses  /A-=l  for 
the  iteration  solutions  of  formulas  (1)  and  (2) .  At  the  same 
time,  calculated  leading  adge  equation  (5)  and  leading  edge 
boundary  conditions  (3)  thus  obtain  the  finally  needed  calcu¬ 
lation  results. 

We  can  see  from  this  that  the  main  calculation  quantities 
are  carried  out  by  the  TSD  equation  and  only  in  the  final  stage 
do  we  use  the  TSDH  equation  and  leading  edge  equation.  There¬ 
fore,  the  calculation  method  in  this  paper  is  much  simpler  than 
directly  solving  the  FVP  equation. 

IV.  Calculation  Results  and  Discussion 

We  used  the  above  mentioned  method  to  calculate  the  pressure 
distribution  of  the  ONERA  M6  wing  under  two  types  of  situations. 
The  leading  edge  sweepback  angle  of  this  wing  is  30°  and  the 
aspect  ratio  and  taper  ratio  are  3.8  and  0.562,  respectively.  Since 
this  wing  has  a  pointed  D  airfoil  section,  the  form  of  its  pressure 
distribution  is  rather  complex  and  under  certain  flow  conditions  it  has 
a  complex  shock  wave  form.  Therefore,  this  wing  has  become  the  standard 
model  for  testing  transonic  calc  methods. 

Figure  1  shows  the  TSDH  solution  for  sectional  pressures  at  two  span 
wise  stations  under  supercritical  conditions  with  no  shock  wave  when 
vj^=  0.70  and  a  =  3°,  and  it  is  very  close  to  the  wind  tunnel  tests  (7) 
and  FVP  solutions  (8>*  8 


Figure  2  gives  the  section 
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pressure  distribution  an  six  spanwise  positions  when  =0.84 
and  a=3.06°  and-  it  is  compaxad_wi4h  the  test  values  [7]  and 
FVP  solution  [Q] .  We  can  see  from  the  figure  that  the  TSDH 
solution  is  very  close  to  the  test  values  and  FVP  solution  and 
relatively  well  reflects  the  complex  shock  wave  form  on  the 
wing  surface.  On  a  section  with  z=80%,  the  TSDH  solution  is 
even  closer  to  the  test  values  as  compared  to  the  FVP  solution. 

In  the  z=95%  area,  the  TSDH  solution  has  a  certain  difference  with 
the  test  values  and  the  main  reason  is  there  is  no  further 
division  of  the  grid  between  z=95%  and  the  wing  tip  and  this 
affects  the  accuracy  of  the  results.  The  figure  always  gives 
the  TSD  solution  with  JU,  =0  and  comparison  with  the  TSDH  solution 
shows  that  the  level  of  congruence  between  the  TSD  solution  and 
test  values  is  relatively  poor  and,  cannot  very  well  reflect  the 
complex  shock  wave  form  on  the  wing.  Figure  3  gives  the  shock 
wave  position  on  the  M6  wing  when  =0.84  and  a=3.06°  and  it 
can  be  seen  that  the  TSDH  solution  is  very  close  to  the  test 
values.  After  z^  80%,’  the  calculated  positive  shock  wave 
position  is  more  towards  the  front  than  the  test  values  and  this 
is  because  the  spanwise  grid  near  the  wingtip  is  very  rarefied. 
Therefore,  in  order  to  attain  accurate  results  it  is  necessary  to 
guarantee  the  density  of  the  grid.  When  calculated  in  a  situation 
wherein  =0.84,  in  order  to  approximately  consider  the 
disturbance  of  the  shock  wave  and  boundary  layer,  we  used  the 
non-conservation  shock  wave  point  difference  scheme  [3], 

Figure  4  further  shows  the  necessity  of  guaranteeing  the  grid 
density.  It  can  be  seen  from  the  figure  that  when  the  number 
of  mesh  points  on  the  chord  line  is  53,  the  wing  surface  then 
has  a  leading  edge  suction  peak  and  at  the  same  time  the  shock 
wave  position  is  clearly  reflected.  Therefore,  when  calculating 
the  transonic  pressure  distribution  of  the  wing  surface,  it  is 
necessary  to  give  sufficient  attention  to  the  grid  density  on 
the  wing  surface  otherwise  this  can  seriously  influence  the 
accuracy  of  the  calculation  results.  Lastly,  Fig.  5  shows  the 


situation  wherein  when  the  TSDH  equation  and  leading  edge 
equation  (5)  are  jointly  solved,  iteration  error  E  changes  with 
iteration  number  n.  It  can  be  seen  from  the  figure  that  low 


relaxation  coefficient  Cu  has  very  good  effects  for  quickening 
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the  convergence  speed. 
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Fig.  1  Comparisons  of  pressure  distributions  and  wind  tunnel 
tests  for  M6  wing  when  =0.70  and  a=3°. 

Key:  (l)-(2)  Solution;  (3)  Wind  tunnel  test. 
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Fig.  2  Comparisons  of  pressure  distributions  and  wind  tunnel 
tests  for  M6  wing  when  =0.84  and  a=3.06°. 

Key:  ( 1 ) — ( 3 )  Solution;  (4)  Wind  tunnel  test. 
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Fig.  3  Comparison  of  shock  wave  positions  of  M6  wing  when 
=0.84  and  a=3.06°. 

Key:  (1)  Wind  tunnel  test;  (2)  TSDH  solution. 
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Fig.  4 


Influence  of  grid  density  on  chordwise  pressure 
distribution  results  (TSDH  solution)  for  M6  wing  when 
M  oc  =0.84  and  a=3.06°. 

Key:  (1)  14  points  on  chord  line;  (2)  27  points  on  chord 
line;  (3)  53  points  on  chord  line. 


V.  Conclusion 


The  TSDH  equation  and  simplified  leading  edge  equation  as 
well  as  its  corresponding  boundary  conditions  in  the  three- 
dimensional  flow  used  in  this  paper  can  relatively  well  attain 
the  transonic  pressure  distribution  of  the  wing  and  can  relatively 
well  reflect  the  complex  shock  wave  form  on  the  wing.  Because 
the  TSD  equation  does  not  consider  the  cross  derivative  of  the 
velocity  potential,  its  calculation  accuracy  is  relatively  poor. 
Because  this  paper  does  not  carry  out  coordinate  conversion,  it 
can  therefore  conveniently  consider  the  complex  configuration  of 
the  wing  as  well  as  the  fence,  leading  edge  sawteeth  and  wingtip 
winglet  etc.  aerodynamic  equipment  on  the  wing  surface;  however, 
the  nonadvantageous  factors  follow  the  densif ication  of  the 
chordwise  division  points  and  it  is  necessary  to  correspondingly 
and  suitably  increase  the  number  of  chordwise  division  points 
otherwise  it  easily  causes  fluctuation  of  the  calculation 
results  of  the  separate  points  near  the  leading  edge. 

This  research  work  received  the  enthusiastic  attention  and 
guidance  of  Professor  Luo  Shijun  and  we  would  like  to  thank 
him  here. 

Appendix  A  Conditions  of  Shock  Waves  Existing  on  an  Infinite 
Swept  Wing 

Reference  [1]  discussed  the  conditions  of  shock  waves 
existing  on  an  infinite  swept  wing.  If  the  sweepback  angle  of 
the  shock  wave  or  leading  edge  of  the  wing  is  _A  ,  then  the 
conditions  of  the  existing  sweepback  shock  waves,  which  is  the 
velocity  component  of  the  shock  wave  normal  line  direction  on 
the  wing  plane,  must  be  greater  than  the  local  velocity  of  sound. 
It  can  be  expressed  as 

a  (Al) 
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In  the  formula,  x  is  the  chordwise  coordinate  of  the  infinite 
swept  wing  (see  Fig.  A-l) . 


Fig.  A-l  Coordinate  direction  of  infinite  swept  wing. 

Because  the  flow  component  does  not  change  along  the  z  axis, 
formula  (Al)  can  be  rewritten  as* 


<J>,<4>S<  (<!>«)«»  CA2) 

in  the  formula,  *  £  ¥  vn  and  (  4>  v) 


corresponds  to  the 


x  w  •  xo  .  x  max 

maximum  value  of  a=0.  From  the  energy  equation,  we  can  obtain 


(^*)  ■»»*  —cos  A 


1  + 


Y  -  1 


(.Vfioos*A)’1 


(A3) 


Assuming  =0 ,  we  can  obtain  different  <p  expressions  for 

y  * 

the  different  velocity  potential  equations.  For  the  FVP  equation 


« - cos*a|i  1  - 7^TCl-(-V/ico3*A)*‘)J 


(A4) 


For  the  TSD  equation 


x,  sec1!  —  A/1 
(  Y  +  1  ).V/i 


(A3) 


For  the  TSDH  equation 


«=*C»/fli-2/l*(.V/i-sec,A)-5JM 
In  the  formula 

A*-(Y  +  1  )A/l+4.Vit5s.V 
S,*(Y+  1  )-V/isec*A 
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(AS) 


For  the  longitudinal  large  disturbance  velocity  potential 
equation  [5] 


/  1  2  (Ml- sec'-\)  _  . 

1 - A.  1 


In  the  formula 


(A7) 


A.- C  Y  + 1) Ml+  (  y 

Figure  A-2  gives  the  relationship  of  the  <p  obtained  from 
formulas  (A4)-(A7)  when  =0.85  which  changes  with  sweepback 

angle  -A .  it  can  be  seen  from  the  figure  that  for  an  infinite 
swept  wing  with  ^50°,  the  TSD  equation  is  unable  to 
obtain  shock  waves,  the  longitudinal  large  disturbance  equation 
causes  -A  to  increase  to  60°  and  the  TSDH  equation  causes  A 
to  increase  to  80°. 


Fig.  A-2  Conditions  of ^existence  of  shock  waves  on  infinite 
swept  wing  <px<  $  ^  ^  for  =0.85. 

Key:  (l)-(2)  Equation;  (3)  Longitudinal  large 
disturbance  equation;  (4)  Equation. 
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A  NUMERICAL  SOLUTION  OF  THE  SHOCK-TURBULENT  INTERACTION  OVER 
A  COMPRESSION  CORNER  FOR  SUPERSONIC  FLOW 

Cao  Qipeng 

Nanjing  Aeronautical  Institute 
Abstract 

This  paper  calculates  the  shock-turbulent  boundary  layer 
interaction  over  a  compression  corner  for  supersonic  flow.  The 
calculations  use  the  Cebeci-Kellet  Box  method;  the  turbulence 
model  uses  the  algebraic  eddy  viscosity  model;  the  pressure 
strength  distribution  uses  the  unified  hypersonic  and  supersonic 
formulas  of  flowing  over  a  wedge;  we  carry  out  iteration 
correction  of  the  shock-turbulent  boundary  layer  disturbance. 

The  calculations  predict  the  wall  pressure  distribution  as  well 
as  the  initial  rise  point  location  of  the  pressure. 

I.  Introduction 

The  compression  corner  for  supersonic  flow  is  a  classical 
example  of  interference  between  the  shock  wave  and  boundary  layer 
The  turbulent  flow  especially  has  even  more  real  significance. 

It  is  related  to  the  operating  efficiency  when  an  aircraft 
deflects  the  operating  surface  and  it  is  also  related  to 
supersonic  inlet  surge.  Thus,  it  has  consistently  been  an 
important  problem  which  has  gained  a  great  deal  of  attention  in 
aerodynamic  design. 

Early,  in  the  1950's,  Drougge  et  al  [1,2, 3, 4]  carried  out 
experimental  research  on  the  problem  of  disturbence  between  the 
shock  wave  and  boundary  layer  produced  by  the  supersonic  flow 
over  a  compression  corner.  The  earliest  theoretical  calculations 
used  the  boundary  layer  momentum  integral  relation  method  and  it 
was  quite  rough.  For  improvement,  in  the  1970's,  they  began 


using  a  succession  of  methods  [5,6]  and  at  the  same  time 
following  the  advances  in  aerodynamic  numerical  calculations, 
two  other  calculation  methods  appeared.  One  type  used  the 
finite  difference  method  to  solve  the  boundary  differential 
equation,  carried  out  iteration  coupling  with  external  inviscid 
flow  [7]  or  carried  out  boundary  layer  inverse  problem  solution 
of  the  disturbance  region  [8].  This  method  separately  solves 
non-separated  and  separated  flows  and  when  compared  with 
another  method  of  directly  numerically  solving  the  Reynold's 
mean  Navier-Stokes  equation  [9,10,11],  the  calculation  accuracy 
is  a  little  poorer  especially  for  separated  flows  where  it  is 
necessary  to  make  large  improvements.  However,  it  has  a  small 
amount  of  calculations,  it  saves  on  computer  time  and  it  is 
even  more  suitable  for  the  present  level  of  China's  computers. 
This  paper  uses  the  first  type  of  method  under  these  assumptions 
yet  uses  the  Cebeci-Keller  Box  method  for  solution  of  the 
equation,  the  further  simplification  of  calculations  and  the 
saving  of  computer  time.  Accuracy  is  relatively  good. 

t 

II.  The  Cebeci-Keller  Box  Solution  Method  For  Binary  Constant 
Compressible  Turbulent  Flow  Boundary  Layer  Equations 

1.  Equations  and  Conversions 

The  set  of  binary  constant  compressible  turbulent  flow 
boundary  layer  equations  are 


Continuous 


-^r(p2)  +  -^r(p5)-0 


(1) 


Momentum 


pu 

Energy 


,  du 

^  Afi  , 

dU 

dp 

,  *  (, 
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*y  V 
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+  p, 
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(2) 


(3) 
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In  the  formula,  -p  u'  v'  =  u  is  called  the  Reynold's  stress 

and  based  on  the  Boussinesq  eddy  viscosity  formula  it  can  be 
written  as: 


,  t,  (i) 

ay 


Key:  (1)  Is  the  eddy  viscosity  coefficient. 


In  the  same  way,  -  P 


—  ir= 


Prt  ay 


Pr^  is  the  turbulent  flow 


Prandtl  number  which  this  paper  takes  as  0.9;  Pr  is  the  molecular 
Prandtl  number  which2this  paper  takes  as  0.72.  H  represents  the 
total  enthalpy,  H=  ~-+h,  h  represents  the  static  enthalpy;  the 
above  each  physical  quantity  indicates  the  mean  number  and 
u',  v*  and  H'  indicate  the  pulsation  values  of  the  velocity  and 
total  enthalpy. 

The  boundary  conditions  are : 

5(*,0)-0,  o(x ,  oi*0,  lim  «(*r  y)-a.(x) 

(1) 

H(*,0)-S-  m(#)m  (° 

Key:  (1)  Or. 

Lower  symbol  w  indicates  the  wall  surface  value;  e  indicates  the 
boundary  layer's  peripheral  value. 

We  carry  out  conversion  of  the  above  set  of  equations  and  let 


From  the  continuity  equation  we  can  introduce  the  flow 


function 


^-(p,n,a^e)1/1/  (6,  n) 


Substituting  in  formulas  (2)  and  (3) ,  after  arrangement  and 
obtaining  conversion  the  equations  are: 

Momentum 

~f,Jk )  (7) 

Energy 

+«//'/')'  +/’./*'  “(/'  -*)  (  8  } 

In  the  formulas 


*  b-e  (  1  +  e«/*)»  c  -pP/P*»»« 

e_»e,/P,  v-jl/p 


The  sign  carrying  a'/)  'indicates  the  derivative  of  rj  . 

The  boundary  condition  conversions  are: 

1  ■  0  i  /;  ■  0 ,  /*0t  *  H»/ H <3s/  (  9  ) 

H  lim  /,  ™  1 1  lira  *  "  1  d) 

1  ■*«  1 

Key:  (1)  Or. 

Formulas  (7)  and  (8)  above  which  solved  the  velocity 
distribution  and  enthalpy  distribution  (or  temperature  distribu¬ 
tion)  in  the  boundary  layer  also  necessitates  the  addition  of 

relational  formulas  for  the  a,,  p  and  £  distributions  in  the 

r  m 

boundary  layer.  Here,  the  viscosity  coefficient  uses  the 
Sutherland  formula 


*-0yfii^**7*s  a)  (10) 

0  -1.4S8X  10*79. 80665-& 

(2) 

2  2 

Key:  (1)  Kg* seconds/m  ;  (2)  Kg* seconds/m  . 


We  use  state  equation  p=  p  RT  to  find  the  density  equation: 

dp  dp  d1! 

Because  of  the  =  ~Ty'=c0’  therefore  the  above 

formula  can  be  changed  to 

P.—pRt 

We  use  the  Cebeci  algebraic  eddy  viscosity  model  to  find  £ 
that  is,  on  the  inner  layer  of  the.  boundary  layer: 


(Oi-P 


au 

*y 


t  “  0  •  ( J  ]  -J-dr]. )  (  1  ■ -  ex  p(- -jfW1  J  0  ^rdT]  )j  /'  (12  a  } 

In  the  formula,  mixing  length  l=ky [1-exp (-y/A) ] 

)"W 


*•“41 


H 

1  —11. 8^— 

p*-. _ 

v,a,  </a. 

r  = 

a?  dx 

,i/» 


On  the  outer  layer  of  the  boundary  layer: 


(e.)„-0.0168  f°°  (a.-  a  )c^y-0.0168v.;?ei',  f°°  (1 


(12  b; 


2.  The  Cebeci-Keller  Box  Solution  Method  [12,13] 


A.  We  introduce  a  new  variable  and  change  formulas  (7)  and 
(8)  into  a  set  of  five  first  order  nonlinear  equations.  Letting 


g-f'  03) 

u>i~g'  (14) 

i  —s'  (15) 

then  formulas  (7)  and  (8)  change  into 

(bw)'  +Plfw+Pi  [~y — (16) 
(et  +  dgui)'  +PJt~  —  t  (17) 


The  boundary  conditions  are: 

TJ  »  0  •  0  y  Qmt  *  Of  $  ** 

0  -T)-,  5-  1,  s  =  1  (1) 


(18) 


Key:  (1)  Or. 


B.  Difference  Equations 


When  using  the  difference  for  solution,  the  grid  layout  is 

.-Ui+if*  «-1.2, . ,JV 

0.  —  0  ,  *1/ —  H/-i  +  hjy  }  —  1 ,  2  , . •  / 

^=■0- 

In  order  to  make  the  near  wall  grid  dense  and  distant  wall 
grid  rare  'ed  so  as  to  suit  the  needs  of  boundary  layer  calcula¬ 
tions  and  not  cause  an  increase  in  the  number  of  direction 
grid  points,  we  used 


We  consider  that  formulas  (13)  ,  (14)  and  (15)  do  not  contain 

x  direction  derivatives  and  therefore  the  Y)  direction  deriva¬ 
tive  uses  the  central  difference  of  the  nodal  point  (  %  n, 
q  j  1/2^  t*ie  ^ n  station  area.  The  function  value  of  the 
nodal  point  uses  the  mean  value  in  the  two  top  and  bottom  points 
and  thus  the  three  difference  equations  are 


(19) 

mi*  «l.  ty* 

A)  (5;  —  <7j-i)  £ 

(20) 

(21) 

.n 


In  the  formulas,  f .  indicates  the  f  value  of  the  q  .  area  in  the 
£  n  station  area  and  the  remaining  are  analogously  deduced. 


Formulas  (16)  and  (17)  possess  q  direction  derivatives  and 
therefore  the  derivatives  are  replaced  with  the  central  difference 
of  the  central  point  (  %,  n_^/2 '  j-1/2^  and  the  function  value 


uses  the  mean  value  of  the  four  points  near  the  central  point, 
that  is,  the  difference  equations  are: 


x  ;*(«?-,  -ir.llK-.'W.ut-fV-li)  3  (23) 


in  the  formulas 


-ty-, 


( -  C(^)?+  (s1)’*1  +  (51)’-,  +  ( 9'-)V-\y  4 
(bw)?1*- ((&»)’+  (W*  V  2 
(ety  »h,-'Uet)r'<'-{ety,:Yt1 


The  remainder  are  analogously  deduced. 


The  boundary  conditions  are: 

/:-  o,  ri-o,  g'~  1 
(1) 

5/*  1 

Key:  (1)  Or. 


(24) 


C.  Numerical  Solution 

If  we  calculate  beginning  from  the  \  -  j|n_^  station,  then 
9j_1/  wj_1,  sj  ^  anc*  ^  are  ^nown  (0  ^  3  ^  J)  and  it  is 

necessary  that  the  f  ” ,  g^,  w^,  s”  and  t?  values  in  the  1,  n 

station  area  be  a  total  of  (5J+5)  and  equations  (19) -(24)  form  a 
set  of  equations  of  these  (5J+5)  unknown  quantities.  In  order 
to  simplify  calculations,  we  broke  them. up  into  two  sets.  We 
used  the  corresponding  boundary  conditions  in  the  first  set  of 
equations  (19) ,  (20)  and  (22)  as  well  as  (24)  to  first  solve  the 
(3J+3)  unknown  quantities  f?,  g"  and  w^;  afterwards,  we  used  the 
remaining  boundary  conditions  in  the  second  set  of  equations 
(21) ,  (23)  as  well  as  (24)  to  solve  the  (2J+2)  unknown  values 

s1?,  and  t1?  and  then  used  the  new  s1?  and  t1?  values  to  replace  the 
3  3  3  3 

coefficient  values  of  the  first  set  of  equations  and  again  solve 
the  new  f?,  g?  and  w?  values.  It  is  repeated  in  this  way  until 

we  satisfy  certain  convergence  conditions. 


For  the  first  set  of  equations,  we  first  separate  and 

classify  the  function  values  of  the  £  and  %,  .  stations  in 

n  n— l 

formula  (22)  into  equal  formulas  on  the  left  and  right  sides  and 
afterwards  use  the  Newtonian  normal  line  to  change  this  set  of 
equations,  that  is: 


We  introduce  iteration  of 


l — 0 ,1,2.  .  .  . 


V  *'  -  -  /  .•  /  V  .  V  / 


The  initial  value  is 

/i5,=  0,  5ot0,*0,  tt-r-wr1 

/;6j=*/r,  gr=9nrl,  i»r=*rl  o  <j</-d 
/r=/7-,  sr  =  i,  <’=u>r 

For  high  order  iteration  we  let 

/r'^fr+b/r 
g^^gf  +  bg?' 
u-r^-wf +6u.’,"> 

Substituting  in  the  first  set  of  equations  and  omitting  the 
second  order  terms  of  £  f ,  5*g  and  <5  w,  we  then  obtain  the 
set  of  linear  equations  (for  convenience  of  calculations,  we 
omit  the  right  upper  corner  (i)  of  each  term): 

bfj-hfi.i — 4f-(65,+6$;-,)«Y  ij  (25) 

6fl,—  65m  — y-(6tu,+  6t»;.i)  -  Y3,  (26) 

(sl);&iu;4-  ($2);&!B/.(t  (s3)/5/;+(44);6//.i+  (s5)/&5/+  (i6)/6<7f.i=“  Y2j  (27) 


In  the  formulas 


VI / — / — /f>  -h  A,p/£i/2 

yzi-gii'i-gr+hiuliU 

Y2/-/?lj:i. -6/ii»/») -h ((-P^.+a.X/ti;)  ' 

+  (CP.).+a.)  (p*)/ii/t+a.(a>j:r.t/;:*/.-/v}/1a»)»1/1)) 

(«i),=A;*6f +-|-C((P1).+o.)/;,>--aJ7:i/1) 

(J2)/  =  -h}'b}i\  +  -i_c((P]).+  «.)//»  -  a 

(J3)/-'|-C((/,1).+a.X,+a,a,v}/1) 

(54) ,«  _1_  C  (  (^,).+ am)u>}t  >,  +  a.u»;:  J/t) 

(*5)y- -((/’.).+ a.)^/0 
(s6)/--C(/>i).+a.)p;» 

^l":I/i--A;‘C(6ui)j-‘-(6w);:l)  +  c-(PI).+  a.)(/u;);:i,.-.(;_(/>s)> 

TO.)(ff!);:!,!-(P!).(P,/p);.l/1-(P1)fcl(pf/p).;;/j 
a»“4«-i/s/ (t«"*4».]) 


In  the  calculations,  the  first  time  we  substitute  the  flow 

parameters  in  the  £  =  .  station  area  into  (  £  )  .  and  ( £'  )  . 

r  n-1  mi  m 

b  is  found  from  the  formula  and  we  let  (  £  ) .=(  £  ).  We  find 

mi  m 

the  yc  value  of  the  intersecting  area  of  the  inner  and  outer 
layers  of  the  boundary  layer.  The  boundary  conditions  remain 
unchanged  in  the  iterations,  that  is,  we  take 

t>/»h ■■  0 ,  bgih  =■  0 ,  bgf  -  0 ,  »-0,  1,  2,  . 

The  set  of  equations  (25) ,  (26)  and  (27)  form  a  block  of  three 

diagonal  matrices  and  by  using  matrix  algebra  we  can  solve  the 

perturbation  quatities  (  6  fi"^  ,  &  g!^  and  w!^).  We  add 

3  D  1 

them  to  the  above  one  time  iteration  value  until  calculation 
convergence  and  the  convergence  data  takes  l<$  w  I  ^  <£  .  The 
lower  symbol  w  indicates  the  wall  surface  value  and  £  is  the 
assigned  small  quantity. 

This  paper  calculates  beginning  from  \  =0,  the  right  side 

of  formula  (16)  is  zero,  the  flow  is  assumed  to  be  a  laminar 
flow  (  ^m=0)  and  we  can  suppose  an  iteration  solution  of  the 

inttial  velocity  distribution  and  enthalpy  distribution  (or 

temperature  distribution) ;  the  next  step  calculates  the  parameters 

of  each  %  station  based  on  the  turbulent  flow  (  £  £0) . 

m 

After  finding  a  preliminary  solution  of  the  set  of  momentum 
equations,  we  can  solve  the  energy  equation  and  following  the 
above  procedure  classify  each  function  value  of  the  if  and 
^  stations  in  the  set  of  equations  in  (23)  to  the  two  sides 

of  the  formula  being  equal  and  write  the  set  of  linear  equations 
as 


(«!!)/*/+  (J12)///.,  +  (^13)5/+  (514)  /5/.,  -  (i?12)/ 


(28) 


In  the  formula 


(sll),= 

(512) /* 

(513) /- 

(514) /- 
(Rl2)i 


■  — |—((^*l).+ a.) - Y~aJV-\ti 

-A;1«f?.t+-|-((P1).+cO/7-i. .  — t 

-a^V-\l 

.-A;,C(«0"*,-(«0":1.)-  2  h-,lt(dgv))Yux-(dgwyi:\,i'y 
-(P,)fc,(/0?:}/S-  2  a^V-VWr-l .i+vJV-ltifV-lt-fi-ifd 


When  the  %  n_1  station  momentum  equation,  energy  equation  as 

well  as  %  n  station  momentum  equation  solutions  are  known,  the 
above  coefficients  (sl2)j  etc.  can  be  found  and  equations  (21) 
and  (28)  and  the  corresponding  boundary  conditions  form  a 
block  of  three  diagonal  matrices  and  we  find  s?  and  tj.  This 
paper  calculates  beginning  from  the  ^  =0  area,  this  area's 
an=0  and  thus  equation  (28)  is  simplified  into: 

(sll)///+(*Z2)/#/.,-(i?12)/  (29) 

In  (sll)  .  and  (sl2)  .,  we  also  omit  the  terms  containing  a  . 

3  1  n 

At  this  time,  under  adiabatic  wall  or  given  t=t  ,  the  first 

w 

block  element  is  singular  and  thus  the  block  element  must 
carry  out  different  processing  and  arrangement  (see  Reference 
[12])  . 


III.  Shock  Wave  and  Boundary  Layer  Disturbance  Corrections 

As  mentioned  previously,  when  there  is  transonic  flow  over 
a  compression  corner,  the  inverse  pressure  gradient  along  the 
flow  direction  caused  by  the  shock  waves  formed  in  the  area  of 
the  compression  corner  causes  the  boundary  layer  in  front  of  the 
angular  point  to  become  thick  which  changes  the  effective  form 
of  the  physical  surface  and  thus  changes  the  pressure  distribu¬ 
tion  along  the  periphery  of  the  boundary  layer.  It  conversely 
also  influences  the  development  of  the  boundary  layer  and  it 
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repeats  thusly  until  establishing  an  equilibrium.  In  this 
way,  for  still  unseparated  flows,  we  first  calculate  the 
inviscid  pressure  distribution  over  a  compression  corner  and 
based  on  this  solve  the  boundary  layer  equation  and  determine 
displacement  thickness  $  *  as  well  as  its  changing  rate  ^ — 


which  changes  with  x.  This  changing  rate  is  the  approximate 
value  of  the  streamline  diflection  angle  induced  from  the 
thickening  of  the  boundary  layer.  This  is  also  the  ratio  of  the 
induced  normal  velocity  and  velocity 
it  is  written  as: 


Wc  on  the  periphery  and 


v 

a. 


•9t* 


db' 

dx 


If  the  included  angle  of  the  original  physical  surface's  form 
and  the  distant  front  overtaking  flow  is  Os,  then  the  included 
angle  of  the  physical  surface's  effective  form  and  distant 
overtaking  flow  after  viscous  disturbance  is 


fc-«-+!3r 


(30) 


Having  &Tf  using  the  Van  Dyke  unified  supersonic-hypersonic 
small  disturbance  theory  [14],  we  can  find  the  new  pressure 
distribution 


P*  1  J.  QJ 

r  y  +  i  / 

fy  +i\ 

1  |  1  ] 

yMi.fi.  ~  y Ml  '  0x 

^  4  ' 

l  4  ) 

1  (Ml- 1)9}  J 

(31) 


Pressure  gradient  parameter  P ^  is 

P  *  du,  x  dpx  d 6r 

*“  a,  dx  “  P^»J  tf0T  dx 


(32) 


The  other  parameters  on  the  periphery  of  the  boundary  layer  are 
all  found  with  p  based  on  the  isentropic  relationship. 


IV.  Calculations  and  Results 


This  paper  carried  out  calculations  with  =2.85, 

Tq=268K,  Re^  /L=6 . 7xl0^/meters ,  adiabatic  wall  surface  and 

compression  corner  a=16°.  The  x  direction  uses  25  grid  points 
and  the  rj  direction  grid  points  number  37  from  the  =0 
area  to  52  grid  points  in  the  final  %  station  area.  For  the 
estimation  of  turbulent  flow  boundary  layer  thickness  %>  ) 

which  enlarges  with  %  used  here  see  Reference  [12] .  The 
^  direction  grids  are  equidistant  and  A  %>  =0.04;  the 

<7  direction  grid  is  the  closest  distance  from  the  wall, 

A  jy  ^=h^=0.01,  m=1.14,  A/)  2=h2  etc*  are  calculated  according 

to  the  above  formula  h.=mh.  . .  The  angular  point  is  located 

3  3  “  4. 

in  the  x=0 . 5  area,  that  is,  between  Nx=13  and  Nx=14,  and  Nx  is 

1 

the  number  of  the  x  direction  station  location.  The  g — tt~  of 

fcyWc 

the  (  Pc/i'Q)  term  in  uses  the  mean  value  of  the  two 

frontcan§  back  angular  points;  when  we  consider  the  existence  of 
the  boundary  layer,  use  of  the  pressure  strength  rise  of  the 
shock  waves  is  not  sudden  but  is  completed  in  any;  multiple 
distance  along  the  wall  surface  equal  to  the  thickness  of  the 
local  boundary  layer.  For  this  reason,  when  beginning  calcula¬ 
tions,  we  complete  the  levelling  of  the  angular  point's  front 
and  back  pressure  rise  in  the  four  x  direction  grids.  In  order 
to  simplify  the  amount  of  calculations,  when  %  ?  0,  we  use  the 
relational  formula  of  the  enthalpy  and  velocity  distribution  when 

P  *1 
r 


Substituting  the  calculation  of  the  energy  equation,  coefficients 
a,b  and  c  are  determined  by  the  adiabatic  wall's  boundary 
conditions.  Calculations  on  the  20,000  times/second  FacoM  230 
computer  only  requires  three  times  of  iteration  to  be  able  to 


obtain  relatively  good  results.  The  C.P.U.  calculation  time 
without  carrying  out  iteration  is  33  seconds  while  the  C.P.U. 
calculation  time  with  carrying  out  iteration  is  one  minute  and 
20  seconds.  The  appended  figures  draw  the  wall  surface  pressure 
distribution,  friction  damping  distribution  as  well  as  the 
velocity  distribution  of  a  station  location  and  they  accord  well 
with  the  wall  surface  pressure  distributions  of  Reference  [15] . 


(3) 

Fig.  1  Comparison  (with  Reference  [15])  of  surface  pressure 
distribution. 

Key:  (1)  Test;  (2)  Calculation  pointy  (3)  Centimeters. 


x  ■  x/L 


Fig.  2  Cf  distribution  along  compression  surface  (the  angular 
point  is  in  the  X=0.50  area,  L=13.8cm). 

Key:  (1)  Test  point. 
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tt/Ur 
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C  6  )  temperature  profile 
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Fig.  3 

Key:  (1)  x=0.52  area  velocity  distribution;  (2)  x=0.52  area 
temperature  distribution;  (3)  Distant  wall;  (4)  Near 
wall. 
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OPTIMUM  DESIGN  OF  AN  AIRCRAFT  WING  STRUCTURE  WITH  MULTIPLE 
CONSTRAINTS  BY  THE  DIRECT-VISION  CRITERIA  METHOD 

Penned  by  Li  Linang 
Shenyang  Aircraft  Corporation 

Abstract 

This  paper  introduces  an  optimized  computer  program  SAFDOP 
which  uses  the  direct-vision  criteria  method  for  a  wing  surface 
structure  with  multiple  constraints.  This  program  can  separately 
carry  out  analysis,  optimization  and  synthetic  optimization 
calculations  of  stress,  displacement,  flutter  and  static 
aeroelasticity .  We  used  the  alternate  constraint  method  for 
synthetic  optimization  calculations.  Aside  from  the  stress 
constraints  using  the  full  stress  method,  each  of  the  other 
constraints  used  the  uniform  derivative  criteria.  For  practical 
purposes,  we  used  the  element  linking  weight  as  an  independent 
design  variable,  in  the- optimum  design  the  automatic  sieve 
selection  design  variable  and  other  methods  were  used  and  we  also 
considered  the  automatic  generation  of  initial  numerical  data. 

We  used  this  program  to  calculate  a  typical  fighter  wing  and  it 
proved  to  be  feasible. 

I.  Introduction 

In  order  to  raise  the  maneuverablity  of  fighters,  it  is 
necessary  to  use  a  thin  wing  with  medium  aspect  ratio  and 
medium  sweepback  angle.  The  structural  design  of  this  type  of 
wing  must  at  the  same  time  consider  the  strength,  flutter, 
static  aeroelasticity  and  displacement  etc.  requirements;  and  to 
raise  the  thrust-weight  ratio  of  the  aircraft,  it  is  also 
necessary  when  satisfying  these  conditions  to  cause  the  designed 
structural  weight  to  be  minimum.  References  [1,2,5]  also  only 
consider  optimization  calculations  under  one  type  of  constraint 


condition.  References  [3,4]  only  performs  optimization  calcula¬ 
tions  under  two  types  of  constraint  conditions.  The  focus  of 
this  paper  is  to  resolve  and  simultaneously  satisfy  the 
optimization  calculations  of  wing  surface  structures  with  various 
types  of  constraint  conditions.  In  this  way,  we  can  reduce  the 
numerical  data  transfer  errors,  save  on  computer  time,  shorten 
the  design  period  and  reach  the  goal  of  the  structure's  minimum 
weight  design. 

II.  Numerical  Data  Generation 

In  order  to  reduce  the  copying  and  transmission  work  of 
initial  numerical  data,  the  SAFDOP  system  considers  the  automatic 
generation  of  intial  numerical  data  and  information.  Its  contents 
include:  (1)  model  node  point  autQmatic  numbering;  (2)  generation 
of  nodal  point  coordinates;  (3)  generation  of  element  information 
included  in  element  linking;  (4)  the  interchange  of  information 
between  the  element  chain  and  element;  (5)  generation  of  the 
quality,  materials,  structural  parameters  and  the  permitted 
stress  etc.  information. 

III.  Analysis  of  the  Structure 

1.  Analysis  of  Stress 

The  analysis  of  stress  uses  the  finite  element  displacement 
method  and  the  displacement  vector  is  obtained  by  the  following 
formula  which  uses  the  variable  band  width  elimination  method. 

(0 

With  the  displacement  of  the  element,  we  find  the  element's 
inner  force  based  on  the  following  formula 

(2> 

We  then  find  the  stress  of  the  element  based  on  the  relation¬ 
ship  of  the  inner  force  and  stress.  Within  this,  b)  is  the 
load  array;  [K]  is  the  structure's  total  stiffness  matrix;  [k]  is 


the  structural  element's  stiffness  matrix;  { s}  is  the 
structure's  inner  force  array;  and  {S}  is  the  structure's 
displacement  array. 

2.  Analysis  of  Flutter 


The  analysis  of  flutter  uses  the  following  equation 


iV 

-co*  2  A/rfft+aJCl  2j  arf9* 

fl-  1  R  *  1 

K-2*(bJhyp 


In  the  formula,  6l>  is  the  flutter  frequency;  q  is  the  generalized 

coordinate;  bQ  is  the  half  root  chord  length;  g  is  the  structure's 

damping  coefficient;  N  is  the  number  of  the  degrees  of  freedom; 

<£  n  is  the  natural  vibraton  frequency;  anR  is  the  generalized 

aerodynamic  force;  M  is  the  generalized  mass;  k  is  the  reduced 
b  CO 

frequency  (k=— - ,  v  is  the  flutter  velocity) ,  footnote  n, 


R=l,2, 


,  N. 


The  sturucture's  intrinsic  frequency  and  form  f  are 

calculated  according  to  the  following  formula  which  uses  the 


iteration  method 


r\  I 

,-<n*CCTTl  M  < 

t  \J 
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In  the  formula,  [C  ]  is  the  flexibility  matrix  and 


is  the  mass  matrix. 


[\1 


[C  ]  is  obtained  from  formula  [1]  and  we  use  the  curved 
surface  spline  function  to  interpose  to  each  mass  point;  after 
obtaining  the  intrinsin  form,  we  use  the  curved  surface  spline 
function  to  interpose  on  12  air  flow  direction  tangent  planes 
(each  tangent  plane  uses  4  points) ,  then  we  use  the  cubic  curve 
to  fit  the  form  of  each  air  flow  direction  tangent  plane  and 


36 


provide  non-constant  aerodynamic  force  for  calculations;  the 
non-constant  aerodynamic  force  is  calculated  by  the  subsonic 
kernel  function  method  and  it  allows  the  use  of  9  control 
points;  we  use  the  v-g  method  to  solve  the  flutter  equation. 

3.  Analysis  of  Static  Aeroelasticity 

Analysis  of  aeroelasticity  uses  the  following  equation 

{<*,>  -  {a,}  +  — Pt»*c C*n  C  A  3  { “,>  (  5  ) 

In  the  formula,  a  is  the  distribution  of  the  initial  attack 

o 

angle;  a,.  is  the  final  attack  angle  distribution  after  consider- 
£  ©Y 

ing  the  aeroelastic  effects;  [C  ]  ' the  corner  flexibility 
matrix;  [A]  is  the  aerodynamic  for.ee  effec  coefficient  matrix. 

Aerodynamic  force  effect  coefficient  matrix  [A]  i.s  calculated 

by  the  Woodward  method  and  we  use  the  curved  surface  spline 

YY 

function  to  interpose  the  C  of  the  structure's  nodal  point  area 

to  each  aerodynamic  force's  divided  characteristic  point  and 

©Y 

then  find  the  air  flow  direction's  slope  and  obtain  C 

For  problems  corresponding  to  different  static  aeroelasticity 
we  give  a  and  after  solving  a^  based  on  requirements  we  find 
the  force  or  moment  of  force  and  lastly  find  ratio  K  of  the 
force  or  moment  of  force  considering  the  aeroelasticity  effects 
and  the  rigid  body  value. 

IV.  The  Structure's  Optimum  Design 

1.  Design  Variables 

In  order  to  avoid  the  changes  of  the  structure's  element 
dimensions  after  optimization  not  satisfying  technical  require¬ 
ments  and  at  the  same  time  to  reduce  the  number  of  design 


variables  so  as  to  make  optimization  calculations  in  medium 
and  small  computers  convenient,  the  SAFDOP  uses  the  element 
linking  weight  as  the  design  variable. 


The  element  linking  is  defined  as  the  joining  of  any  inter¬ 
connected  similar  type  elements.  One  element  linking  includes 
any  similar  type  elements  and  the  dimensions  of  these  elements 
are  obtained  from  two  control  dimension  linear  interpolations 
of  the  element  linking. 


2.  Optimum  Design  of  Stress  Constraints 


We  use  the  full  stress  method  for  stress  constraint  optimiza¬ 
tion.  Under  given  structural  layout  and  materials  conditions  as 
well  as  many  types  of  load  effects.,  each  element  at  the  least 
has  the  operating  stress  reach  an  allowable  value  under  a  type 
of  load  effect  or  reach  limiting  dimensions.  The  expression  of 
its  constraint  conditions  is: 


9 1 


Pi*  1 


(1) 


Key:  (1)  Or. 


(6) 

(7) 


In  the  formulas,  G.  is  the  maximum  operating  stress  of  the 

i  number  elements  under  many  types  of  load  effects;  (  G  is  the 

permissable  stress  of  the  i  number  elements,  its  positive  and 

negative  signs  are  determined  by  the  positive  and  negative  signs 

of  G  .  ,  it  can  be  a  constant  value  and  it  can  also  be  the 

imax 

buckling  damage  stress  which  changes  with  the  design  variables 

in  the  iteration  process;  £  is  the  permissable  error  and  is 

predetermined  by  the  designer;  x.,  x.  and  x.  .  are  separately 
c  1  3  i  imax  lmin  r  1 

the  element’s  calculated  dimensions  and  maximum  and  minimum 
dimensions  corresponding  to  the  i  number  design  variables.  In 


order  to  attain  the  above  goals,  it  is  necessary  to  carry  out 
iteration.  The  iteration  formula  is 


u'i'i = f  °K'  -~-Y 


In  the  formula,  VL  is  the  weight  of  the  i  number  elements  (or 
element  linking) ,  K  indicates  the  number  of  iterations,  R  is  the 
superrelaxation  coefficient.  It  is  used  to  regulate  the 
correction  quantity,  it  has  a  relatively  large  effect  on  the 
convergence  speed  and  it  is  mainly  determined  by  the  structure's 
form  and  constraint  conditions. 

The  definition  of  the  element  linking's  full  stress  is  that 
in  each  element  linking,  at  least  two  points  of  stress  not  on 
the  same  tangent  plane  reach  the  constraint  condition  require¬ 
ments  of  full  stress  and  each  of  the  other  points  do  not  allow 
excessive  stress. 


3.  Single  Constraint  Optimization  Calculations  of  the 
Displacement,  Flutter  and  Static  Aeroelasticity 


Single  constraint  optimization  of  displacement,  flutter  and 
static  aeroelasticity  all  use  uniform  derivative  criteria.  We 
take  the  design  variable  to  be  W,  the  structural  weight  of  the 
wing  surface  is  objective  function  F,  g  is  the  constraint 
function  and  then  the  Lagrangian  function  is 

Q(m-F(H')  +  \gOV)  (9) 

The  limiting  value  condition  is 


_  aF(W) 
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As  a  result 


V 


and  because 

dF(lV)  _  , 
ali'i  "  1 

therefore 

' T-**  U)  <>•) 

Key:  (1)  Constant. 

% 

V. 

We  can  know  from  the  above  formulas  that  for  situations 
with  single  constraints,  when  the  weight  partial  derivatives  of  <* 

the  constraint  function  for  each  element  (or  element  linking)  are 
all  constants,  the  structural  weight  reaches  minimum  value. 

Its  constraint  expressions  are' separately : 

For  displacement  constraints 

(n> 

*  (1> 

a  1  >0  (12) 

«*f  *  i 

Key:  (1)  Or. 

In  the  formulas,  &  is  the  generalized  displacement  under 
aerodynamic  load  effects.  It  can  be  the  torsional  angle  or 
bending  angle  of  a  certain  tangent  plane  on  the  wing  or  the 
deflection  of  a  certain  nodal  point  or  the  turning  angle  of  a 
certain  spar;  (  <5  )  is  its  corresponding  permissable  value. 

These  displacement  constraints  can  be  on  nodal  points  as  well 
as  not  on  nodal  points.  The  displacement  not  on  nodal  points 
is  obtained  based  on  the  linear  interpolation  of  the  two 
adjacent  nodal  point  displacements;  x.  and  x.  .  are  separately 
the  element's  maximum  and  minimum  limiting  dimensions  correspond¬ 
ing  to  the  i  number  design  variables. 
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For  the  flutter  constraints 


9  i 


x‘  0 


(1) 
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(13) 

(14) 


Key:  (1)  Or. 


In  the  formulas,  is  the  claculated  flutter  speed;  (v^)  is  the 
flutter  speed  needed  to  be  reached. 


For  the  static  aeroelastic  constraints 


'■“I1  «)l<e 

(1) 

(15) 

a,  -  1  — 0 
Xi 

•*» 

(16) 

Key : 

(1)  Or. 

In  the  formulas,  K  is  the  ratio  of  the  calculated  flexible  wing 
surface  force  (moment  of  force)  and  the  stiff  wing  surface  value 
(K)  is  the  K  value  which  needs  to  be  reached. 

When  separately  satisfying  these  constraint  conditions,  in 
order  to  reach  the  goal  of  uniform  derivatives  needed  for 
carrying  out  iteration,  the  iteration  formulas  is  similar  to 
iteration  formula  (8)  for  full  stress  design,  that  is 


(i) 

Key:  (1)  Objective. 


(17) 


In  the  formula,  (  3  g/  <?  W)  ^  ject^ve  ^"s  ca^ec^  the  objective 

derivative  and  it  has  the  same  value  for  all  of  the  design 
variables.  In  reality,  it  controls  the  parameters  of  the  design 


variable's  changing  step  length  and  based  on  each  step  we 
anticipate  obtaining  constraint  function  increment  g  and 
determine  the  objective  derivative  values. 


For  the  displacement  constraints 

S  -HOJ  +(^)o  J  )  <l8> 

*  1  (1) 

Key:  (1)  Objective. 

For  the  flutter  and  static  aeroelastic  constraints 


(19) 


In  order  to  control  the  step  length,  we  take  g=  , 

n  is  the  positive  integer  given  by  the  designer,  R  always  takes 
1/2  for  the  flutter  and  displacement  constraints  and  1  for  the 
static  aeroelastic  constraints. 

For  the  flutter  and  static  aeroelastic  constraints. 

The  calculation  expression  of  the  objective  derivative  is 
jointly  obtained  from  formulas  (17)  and  (19) : 


(20) 


Key:  (1)  Objective. 


In  the  formula,  n  is  the  number  of  design  variables. 


For  the  displacement  constraints,  we  can  obtain  the 
corresponding  objective  derivative  expression  in  the  same  way 
from  formulas  (17)  and  (18) . 

In  fact,  when  we  only  have  the  main  effect  of  the  elements 
on  the  constraint  function,  the  constraint  function  can  gradually 
become  uniform  in  the  iteration  process  for  its  derivative. 

These  elements  are  called  "critical  elements".  However,  when 
the  constraint  effects  are  not  great  or  have  an  inverse  effect 
(i.e.  when  the  element  has  a  negative  derivative  and  there  is 
an  increase  in  the  structural  weight,  the  constraint  function 
value  is  actually  reduced)  or  for  elements  wherein  the  constraints 
have  major  effects  and  the  dimensions  exceed  the  maximum 
limiting  value,  there  derivatives  will  tend  to  be  not  uniform. 

The  dimensions  of  these  elements  are  determined  by  the  limiting 
geometric  dimensions  and  these  elements  are  called  "non-critical 
elements" . 

In  order  to  quicken,  convergence,  when  saving  on  computer  time, 
during  iteration,  we  use  the  automatic  sieve  selection  design 
variable  method  which  is  the  automatic  elimination  during  each 
iteration  cycle  of  design  variables  which  belong  to  negative 
derivatives  or  have  already  reached  the  smallest  dimensions  and 
the  derivative  is  smaller  than  the  objective  derivative. 

4.  Calculation  of  Synthetic  Optimization 

(1)  Constraint  Conditions 


(2)  Optimization  Method 


The  optimization  calculation  when  the  wing  surface's 
structure  simultaneously  satisfies  the  requirements  of  the  above 
two  types  or  more  constraints  (aside  from  the  limiting  geometric 
constraints) ,  it  is  called  synthetic  optimization  calculations 
or  multiple  constraints*  optimization  calculations.  The 
synthetic  optimization  method  uses  the  constraint  rotational 
method  which  transforms  the  multiple  constraint  problem  into  the 
processing  of  a  series  of  single  constraints.  We  then  use  the 
selection  of  the  major  constraint  to  reach  the  goal  of  reducing 
the  number  of  constraints.  The  specific  method  is  on  the  basis 
of  the  full  stress  design,  we  check  whether  or  not  the  constraints 
are  satisfied  and  if  they  are  this  is  the  state  of  synthetic 
optimum  design.  If  they  are  not  satisfied,  then  we  distinguish 
by  means  of  the  major  constraints,  select  the  most  "dangerous" 
constraints  and  optimize  them.  Afterwards,  we  check  whether  or 
not  the  remaining  constraints  are  satisfied  and  if  they  are  we 
transfer  into  the  full  stress  calculation,  otherwise  we  carry 
out  major  constraint  selection  and  repeat  the  above  process  until 
the  constraint  conditions  are  all  satisfied  and  continue  until 
the  difference  of  the  weights  calculated  from  the  full  stress  of 


the  two  cycles  reaches  the  predicted  value.  That  is,  the 
structure  reaches  the  state  of  synthetic  optimum  design.  At 
this  time,  the  elements  in  the  structure  are  divided  into 
different  types  of  "critical  elements".  For  example,  the 
"stress  critical  elements",  "displacement  critical  elements", 
"flutter  critical  elements",  "aileron  critical  elements"  and 
minimum  or  maximum  geometric  dimension  elements".  All  of  these 
elements  separately  attain  the  requirements  of  each  optimization 
criterion.  That  is,  the  "stress  critical  elements"  attain  the 
requirements  of  full  stress  calculations  and  the  "displacement", 
"flutter"  and  "aileron  efficiency"  critical  elements  all 
satisfy  the  requirements  of  uniform  derivative  criteria. 
Naturally,  when  the  wing  surface  structure  is  in  a  synthetic 
optimum  state,  these  "critical  elements"  do  not  necessarily 
exist  simultaneously  but  the  number  of  critical  elements  which 
appear  is  related  to  the  constraint  characteristics.  However, 
it  is  necessary  to  point  out  that  when  we  enter  in  certain 
constraint  optimization,  the  "critical  elements"  controlled  by 
other  constraints  can  only  increase  the  dimensions  and  only  the 
dimensions  of  its  own  "critical  element"  can  change  up  and  down; 
the  "critical  elements"  can  be  transformed  during  the  optimiza¬ 
tion  process.  For  example,  when  the  "stress  critical  element" 
enters  into  flutter  constraint  optimization,  this  increases  the 
dimensions  and  changes  it  into  a  "flutter  critical  element". 

See  Fig.  1  for  a  flow  chart  of  synthetic  optimization. 


Fig.  1  Flow  chart  of  synthetic  optimization. 

Key:  (1)  Original  data  generation;  (2)  Full  stress 

design;  (3)  Convergence  discrimination;  (4)  Is; 

(5)  Printed  results;  (6)  Is  not;  (7)  Main  constraint 
discrimination;  (8)  Main  constraint  optimization; 

(9)  Is;  (10)  Is  not;  (11)  Are  the  remaining 
constraint  conditions  satisfied  or  not? 

V.  Examples  of  Utilization 

We  will  here  only  give  one  example  to  show  the  use  of  this 
program  to  calculate  many  types  of  wing  surface  structures.  The 
calculation  model  is  as  shown  in  Fig.  2.  For  thin  wings  with 
medium  aspect  ratios  and  medium  sweepback  angles ,  the  root  and 
fuselage  linking  uses  four  fixed  points  and  one  hinge  point  and 
the  form  is  a  multiwall  structure.  Aside  from  the  three  bent 
beam  roots  being  steel  pieces,  the  other  materials  are  all 
aluminum  alloy.  We  used  the  variable  section  isometric  force 
bar  element  and  quadrilateral  cut  plate  element  to  simulate  the 
spar,  rib  flange  and  skin;  we  used  a  symmmetric  cut  trapezoidal 


plate  element  to  simulate  the  spar  and  composite  plate  of  the 
ribs;  we  used  an  equal  section  bar  element  to  simulate  the 
elastic  fixed  connecting  of  the  wing  and  fuselage;  we  used  the 
"fictitious  bar"  with  relatively  large  stiffness  to  simulate 
the  front  and  rear  end  structures.  We  used  symmetrical  condi¬ 
tions  and  the  half  structure  for  calculations.  The.  calculation 
model  had  a  total  of  106  nodal  points  and  318  degrees  of 
freedom;  in  the  flutter  calculations,  we  used  300  mass  points 
and  in  the  aileron  efficiency  calculations  there  were  44 
aerodynamic  blocks. 

We  separately  carried  out  optimization  calculations  of  the 
element  linking's  full  stress  flutter,  aileron  efficiency, 
wingtip  turning  angle  and  their  synthesis.  The  wing  structure's 
weight  is  given  that  the  element  linking's  full  stress  state 
was  451  kilograms,  the  flutter  speed  was  485  meters/second,  the 
aileron  efficiency  was  0.088  and  the  wingtip  turning  angle  was 
a  0.0033  radian.  The  design  required  that  the  flutter  speed  be 
550  meters/second,  the  aileron  efficiency  K  value  be  0.3  and  the 
wingtip  turning  angle  be  0.03.  The  calculation  results  are 
given  in  Figs.  4-7. 


Fig.  2  Model  of  wing  structure. 

Key:  (1)  Shows  the  flexibility  coefficient  point; 
(2)  Fictitious  bar. 


Fig.  3  Typical  model  of  skin  element  linking. 
Key:  (1)  Aileron. 


Fig.  4  Iterative  process  of  flutter  constraints. 

Key:  (1)  Meters/second;  (2)  Required  velocity  theoretical 
increment;  (3)  Kilograms. 


Fig. 


Iterative  process  of  aileron  efficiency  constraints. 
Key:  (1)  Required  theoretical  K  value;  (2)  Kilograms. 


Fig.  6  Distribution  of  skin  thickness  for  full  stress  design 


Fig.  7  Distribution  of  skin  thickness  for  synthetic  optimum 
state. 

In  the  figure:  %  is  the  thickness  increment  of  the 
critical  element  of  the  skin  thickness  relative  full 
stress  design;  Q  is  the  aileron  efficiency  critical 
element;  (|)  is  the  stress  critical  element;  Q)  is  the 
minimum  dimension  element. 

• 

The  simple  constraint  calculation  results  show  (see  Figs.  4 
and  5)  that  it  is  only  necessary  to  have  3-4  iteration  cycles 
to  converge  to  the  optimum x  <esign  state;  the  flutter  speed  rose 
13.6%,  the  weight  increased  3.75%;  the  aileron  efficiency  rose 
275%  and  the  weight  increased  14%.  Synthetic  optimization 
results  show  (see  Fig.  7)  that  it  is  only  necessary  to  use 
two  iteration  cycles  to  reach  the  synthetic  optimization  state, 
the  main  constraint  is  the  aileron  efficiency  constraint  and 
aside  from  the  root  being  the  "minimum  dimension  element”  and 
"stress  critical  element",  all  of  the  others  are  "aileron 
efficiency  critical  elements".  The  trailing  edge  skin  of  the 
aileron  area  and  its  inner  side  (5,9  element  linking)  has  the 
most  noteable  effects  on  aileron  efficiency.  The  thickness  of 
the  "stress  critical  elements"  6,13  and  14  are  all  reduced  as 
compared  with  the  thickness  during  the  first  full  stress.  This 
is  due  to  the  increases  of  the  thicknesses  of  the  other  element 


linkings  which  causes  rematching  of  the  inner  force.  Synthetic 
optimization  results  also  show  that  on  the  basis  of  full  stress, 
we  must  cause  the  K  value  to  rise  from  0.088  to  0.3  (increase  of 
240%)  which  requires  increasing  the  weight  69.5  kilograms 
(increase  of  15.4%);  at  the  same  time,  the  flutter  speed  .rises 
from  485  meters/second  to  592  raeters/second  (increase  of  22%). 

We  can  see  from  this  that  for  the  wing  surface  structural  design 
of  modern  fighters,  it  is  necessary  to  simultaneously  consider 
the  requirements  of  aeroelasticity  and  static  strength. 

VI.  Concluding  Remarks 

The  SAFDOP  system  provides  a  feasible  method  for  synthetic 
optimization  calculations  of  the  wing  surface  structure.  After 
several  iteration  cycles,  we  can  reach  the  goal  of  an  approximate 
minimum  weight  design,  yet  it  is  also  necessary  to  make  further 
developments  and  improvements  so  as  to  adapt  to  the  needs  of  the 
synthetic  optimum  design  of  composite  material  wing  surface 
structures. 

The  SAFDOP  system  computer  program  of  optimum  design  of  a 
wing  surface  structure  with  multiple  constraints  was  jointly 
completed  by  Lin  Liang,  Guo  Feng,  Wang  Shihao,  Zhang  Yining, 

Sun  Qishan,  Ba  Hengyu  and  Yao  Guijun.  We  also  received  the 
guidance  of  high  level  engineer  comrade  Guan  De. 
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THEORETICAL  SPECIFIC  IMPFLSE  FOR  TWO-PHASE  FLOW  IN  NOZZLE  OF 
ROCKET-RAMJET  ENGINE-  A  SIMPLE  ESTIMATION  METHOD 
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Beijing  Research  Institute,  China  Precision  Machinery  Company 
Abstract 


This  paper  derives  a  theoretical  specific  impulse  formula 
for  two-phase  flow  in  nozzle  rocket-ramjet  engine  on  the  basis 
of  similarity  to  that  used  for  rocket  engines.  Theoretically, 
it  is  proved  that  the  properties  of  two-phase  flow  in  nozzle 
for  these  two  types  of  engines  can  be  analyzed  with  a  common 
combined  parameter  which  is  equivalent  to  the  exhaust  velocity. 
We  introduce  a  simple  estimation  method  for  specific  impulse 
losses  of  two-phase  flow  with  variable  lag.  We  propose  an 
approximate  solution  to  deal  with  the  problem  of  singularity. 


:  :  gram  radius 

P 

’  :  temperature 

i  :  velocity 

r  :  specific  heat  ratio 

p  :  density 

r)  :  two-phase  flow  rate  ratio 

i/,  :  gas  viscosity  coefficient  of  combustion  chamber 
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relative  velocity  lag  quality 
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:  relative  temperature  lag  quantity 


:  residual  air  coefficient 


:  theoretical  air-fuel  ratio 


Lower  Symbols 

o  :  equilibrium  quantity  without  lag 

g  :  gas  quantity 

p  :  grain  quantity 

M  :  two-phase  mixture  quantity 

c  :  combustion  chamber  quantity 

s  :  particulate  matter  quantity 

i  :  nozzle  inlet  section  quantity 

e  :  nozzle  outlet  section  quantity 

X.  Use  of  the  Theoretical  Specific  Impulse  Formula  For  Two- 
Phase  Flow  in  Nozzle 

Based  on  the  two-phase  flow  theory  of  solid  rocket  engines 


[1-5],  under  design  conditions,  the  two-phase  flow  theoretical 
specific  impulse  of  the  rocket  engine  nozzle  can  generally  be 
expressed  as  follows: 


1  +  V». 

i 


Can  we  use  a  similar  method  to  calculate  the  theoretical 
specific  impulse  of  a  rocket-ramjet  engine  nozzle?  Truly,  in 
the  rocket-ramjet  engine,  due  to  the  entrance  of  air,  the 
combustion  process  of  its  fuel  and  the  rocket  engine  have  marked 
differences  and  each  type  of  rocket-ramjet  engine  also  have 
differences  because  of  the  air  permeating  into  the  after 
combustion.  Thus,  as  regards  the  flow  of  the  gas  flow  in  the 
nozzle,  because  we  can  always  view  the  gas  mixture  of  the  nozzle 
inlet  area  as  being  formed  by  two-^phase  material,  we  should  also 
be  able  to  use  a  method  similar  to  that  for  rocket  engines  to 
study  the  problem  of  its  two-phase  nozzle  flow.  In  order  to  not 
lose  universality,  we  selected  the  after  combustion  and  completely 
combined  rocket-ramjet  engines  as  the  objects  of  research.  Based 
on  the  law  of  the  conservation  of  momentum,  after  suitable 
conversion,  we  can  indicate  the  fuel  equivalent  theoretical 
specific  impulses  of  these  two  types  of  engine  nozzles  in  the 
following  forms: 

After  combustion: 


Completely  combined: 


-H( 


l  +  *  +- 


*4* 

1  +TloM 


t-na.  j 


It  can  be  known  from  the  above  formulas  that  even  though  the 
operating  principles  of  the  rocket  and  rocket-ramjet  engines  are 
different  and  the  combined  formulas  of  each  type  of  rocket- 
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ramjet  engine  are  different,  yet  the  rules  of  the  effects  of 
the  grains  on  them  are  very  similar  as  both  are  shown  by  a 
common  combined  paramter  u^.  Therefore,  we  can  use  the 

research  on  this  combined  parameter  to  reveal  the  common  effect 
laws  of  the  two-phase  flow  effect  on  the  performance  of  each 

1+/)0^u 

type  of  engine.  —  u^  can  then  be  viewed  as  the  equivalent 

outlet  exhaust  velocity  of  the  nozzle  when  there  is  two-phase 
flow  and  thus  the  problem  of  two-phase  nozzle  flow  specific 
impulse  loss  actually  becomes  the  problem  of  the  two-phase 
nozzle's  equivalent  outlet  exhaust  velocity.  Its  general 
formula  can  be  expressed  as  follows: 


r  /  v  - 1  w 

1  +  V».  „  .  _L±2ag«_  gr  T  i  ^ — 1 

l+i|t  l+n,  L  \  S  ). 


In  the  formula 


(G;<+nl)corC,)/(  i  +IJ0ci)i) 
v  =fl,V,/C  1  +(Br-  1 
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6>;.  Cd  ^  and  are  separately  the  differentials  of  *C>ui  o>T  and 

T  for  %  . 
g  0 

II.  Solution  of  the  Equivalent  Exhaust  Velocity 


Aside  from  the  simple  two  phase  flow  without  lag  and  with 
constant  lag  of  u>  ^  in  most  situations,  because  of  the 

CO  and  t£>  ^  quantities  of  the  changes  produced  along  the 
nozzle's  axial  and  radial  directions,  there  is  no  way  of  using  a 
simple  method  for  direct  solution.  There  are  two  main  approxi¬ 
mate  solutions  for  the  variable  lag  two-phase  flow  process:  one 
is  the  analytical  method  and  the  other  is  the  numerical  value 
method  [6-8].  During  rough  calculations,  we  can  use  the 


u 
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analytical  method  and  although  the  precision  is  a  bit  poor  yet 
it  can  allow  engineers  and  technicians  to  directly  view  the 
dependent  relationship  between  the  various  initial  parameters 
which  they  selected  themselves  and  the  calculation  results. 

They  can  thus  clearly  know  how  the  selection  of  parameters  can 
causes  the  calculation  results  to  face  their  own  hoped  for 
direction  changes.  When  we  use  the  analytical  method  to  find 
the  equivalent  exhaust  velocity,  its  formula  can  be  expressed  as 


1  «  -fit  ^ 

/(C)  1 

1-M,  *  ‘l1  •  l+’l. 

In  the  formula 


1  +-Ii-K(v,-1)A/J 

_ lA _ 


Under  Stokes  flow  conditions,  because  f ,=f  =1,  we  can  then 

d  n 


further  simplify  it  into: 


.  ft*  l+K(v„-lWo  C'  ,4f 

c, - w,  irdL 


When  compared  with  the  two-phase  flow  without  lag,  the 
specific  impulse  loss  created  by  the  lag  is  expressed  by  a 
relative  quantity: 


1  +T]0  V|  M\t 


The  above  formula  shows  that  the  specific  impulse  loss 
created  by  grain  lag  is  not  only  related  to  the  basic  character¬ 
istic  parameters  of  q’  ^  o  '  V0  etc*  tw°~phase  combustion 
products,  but  is  also  related  to  the  size  and  change  law  of  the 
nozzle  form’s  characteristic  parameter 


Thereofre,  under  two-phase  flow,  this  causes  the  other  conditions 
to  be  the  same  and  only  because  the  nozzle  forms  are  different 
this  can  in  the  same  way  cause  the  same  engine  to  have  different 
specific  impulse  losses. 

III.  Nozzle  Form  Coefficient 


In  order  to  compare  the  effect  of  the  nozzle  forms  on  the 
two-phase  flow  characteristics,  we  can  use  a  nozzle  form 
coefficient  f  to  express  the  ratio  of  a  common  form  nozzle's 
specific  impulse  loss  and  the  specific  impulse  loss  of  the 
minimum  loss  nozzle: 


/. 


A7*  /  n. 

*  A/,.,  \1+X 


\  /f  «,  \  JJU 

vJ3/'5,y/  vj M\J  mi*  J  tt.). 


(8) 


When  the  operating  conditions  of  the  powder  charge  composi¬ 
tion  and  the  engine  are  the  same,  the  above  formula  can  be 
simplified  into: 

(8a) 


In  the  formula 


+  *(v,-  1  )M\ 


ii/t  -i 
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In  this  way,  the  effects  of  the  nozzle  form  can  be  made  to 

stand  alone  from  the  total  effeciency  of  the  nozzle  and  be  a 

single  factor  for  carrying  out  deeper  research.  At  the  same 

time,  we  also  provide  a  new  method  for  experimental  research  on 

the  nozzle  form  characteristics  and  comparing  the  differences 

between  the  theoretically  calculated  results  and  the  test 

results.  Calculations  show  that  the  value  of  common  nozzle  form 

coefficient  f  is  approximately  between  1.5  and  2.5. 
z 
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IV.  An  Approximate  Solution  of  the  Singularity 


Because  nozzle  form  parameter  £  '  is  a  quantity  of  sinularity 
existing  in  the  throat,  there  is  generally  no  way  of  directly 
finding  its  thraot  value  and  thus  there  is  also  no  way  of 
directly  finding  the  value  of  arbitrary  function  integral  term 
J  (  £  ) .  In  the  references,  they  proposed  an  approximate  solu¬ 
tion  of  the  singularity  problem  of  several  types  of  non-isentropic 
flows  and  thus  as  regards  the  two-phase  nozzle  process,  some  of  ‘ 
these  solutions  are  not  very  suitable  while  some  have  relatively 
low  precision.  In  order  to  attain  simplicity  as  well  as 
sufficient  precision,  here  we  will  propose  a  constant  slope 
method.  Its  formula  is: 


r,  i  r.r 
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In  the  formula 


%  .  and 
3  im 


^  is  the  relative  pressure  value  of  the  nozzle  throat  area; 
£  eQ  are  separately  the  relative  pressure  values 


of  the  upstream  and  downstream  joining  point  area  close  to  the 
nozzle  throat. 


Institute  Vice-President  Liang  Shouban  went  over  this  paper 
and  Assistant  Head  Engineer  Bao  Kerning  proposed  valuable  views 
for  guidance  and  we  would  like  to  express  out  sincere  thanks  to 
them  here. 
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